% MT fwd & inverse computaton - keller & kaufman 
% latest date 04.01.2002
% inverse problems ----(???)
% fwd nlayer bug fixed 02.1.2002

% t = period (input - HZ)
% rho = layer resisitivity ohm-m
% h = layer thickness (IN KM)
% j9XY;
% added to check the s impedance proposed by Sri Niwas et al
% 26.2.5


function[z1] = s_imp(rho,h,t),

h = h/1000;
n = length(t);
nlayer = length(rho);

if nlayer ~= length(h)+1,
   fprintf('Nlayers & Nresitivities do not match\n exit\n');
   return;
end;
% declarations

      c= 1+ sqrt(-1);
      pi4=pi/4;
      cte=2*pi4*4/sqrt(10);
      amu=1.2566e-6;
      term=1.e-6/(8*pi4*amu);

for it = 1:n,
      
%formulae for id mt

	for i = 1:nlayer,
		b(i)=cte*sqrt(rho(i)/t(it));
		s(i)=c*b(i);
	end;
	z=s(nlayer);


	for i = 1:(nlayer-1),    
		ni=nlayer-i;
        		a=(z-s(ni))/(z+s(ni))*exp(-c*2*h(ni)*b(ni)/rho(ni));
		z=s(ni)*(1+a)/(1-a);
 	end;
      	z1(it)=z;
end;



